Spin and orbital magnetic response in metals: susceptibility and NMR shifts 



Mayeul d'Avezac^, Nicola Marzari^, and Francesco Mauri^ 
^Institut de Mineralogie et Physique des Milieux Condense, 
case 115, 4 plo.ce Jussieu, 75252, Paris cedex 05, France and 
^Department of Materials Science and Engineering, 
Massachusetts Institute of Technology, Cambridge, Massachusetts 02139-4307 
(Dated: February 1, 2008) 

A DFT-based method is presented which allows the computation of all-electron NMR shifts of 
metallic compounds with periodic boundary conditions. NMR shifts in metals measure two compet- 
ing physical phenomena. Electrons interact with the applied magnetic field (i) as magnetic dipoles 
(or spins), resulting in the Knight shift, (ii) as moving electric charges, resulting in the chemical (or 
orbital) shift. The latter is treated through an extension to metals of the Gauge Invariant Projector 
Augment Wave(GIPAW) developed for insulators. The former is modeled as the hyperfine inter- 
action between the electronic spin polarization and the nuclear dipoles. NMR shifts are obtained 
with respect to the computed shieldings of reference compounds, yielding fully ab-initio quantities 
which are directly comparable to experiment. The method is validated by comparing the magnetic 
susceptibility of interacting and non-interacting homogeneous gas with known analytical results, 
and by comparing the computed NMR shifts of simple metals with experiment. 

PACS numbers: 71.45.Gm, 76.60.Cq, 71.15.-m 



I. INTRODUCTION 

Nuclear Magnetic Resonance (NMR) is a widely used 
and powerful technique for structural determination, 
both in chemistry and in solid-state physics^. It also 
yields valuable information on the electronic struc- 
ture of solids. For instance, NMR was instrumental 
in determining the dx2_y2 pairing of high-temperature 
superconductors^. Empirical rules have been determined 
which relate NMR quantities to physical and chemical 
properties. Unfortunately, such rules can become inac- 
curate when subtle quantum effects are involved. In this 
work, we provide a method for computing NMR shifts 
from first-principles in metallic systems with periodic 
boundary conditions. 

Recent advances have made possible the computation 
of NMR shifts in molecules^ and insulating solids with 
periodic boundary conditions^i^, leading to a better in- 
terpretation of experimental data in systems as diverse 
as zeolite^ or vitreous Boron oxides^. 

At present, to the best of the authors' knowledge, there 
is no complete ab-initio theory of NMR shifts in metal- 
lic systems. Indeed, NMR shifts in metals result from 
two different physical phenomenon. The electronic struc- 
ture can react to the external magnetic field (i) as a dis- 
tribution of magnetic spins, giving rise to the Knight 
shift, (ii) as a distribution of electronic charges, with 
the NMR orbital shift as a result. In most metallic sys- 
tems, the NMR shift is dominated by the Knight shift 
contribution, sometimes by as much as two orders of 
magnitude. As such, it has been the subject of many 
theoretical studies^i^. On the other hand, the devel- 
opment of methods capable of computing orbital shifts 
in metallic compounds has been lagging behind. Yet, 
experiments do not distinguish between the shifts aris- 
ing from these two phenomena. Furthermore, experi- 



mental shifts are given with respect to some insulating 
reference-compound. As such, theoretical calculations 
must include both orbital and Knight shifts in the ma- 
terial of interest and a reference-compound before being 
compared to experiment. The Knight shift is related to 
the density of s-states at the Fermi level. As such, there 
are a number of systems for which the Knight and or- 
bital contributions to NMR shifts and to the magnetic 
susceptibility are of similar magnitude. These systems 
include semi-metals such as graphcne, graphit o^'^i^^ , in- 
tercalated graphitei^ii^, and nanotubesi^, metals with 
strong d-character such as Platinum catalystsi^ii^, or or- 
ganic compounds adsorbed upon metallic catalyst o^^i^^ . 

The aim of the method presented here, is to provide 
a unified first-principles framework to compute both or- 
bital and Knight shifts in metallic systems with periodic 
boundary conditions. The setting for the method is den- 
sity functional theory (DFT) as implemented in plane- 
wave, pseudo-potential codes. The projector-augmented 
wave (PAW) approacbi^ allows us to obtain accurate re- 
sults from pseudo-potential quantities. The problem of 
gauge invariancc in periodic pseudo-potential systems is 
treated using the gauge- invariant projector augmented 
wave (GIPAW) approach of Ref. |j. Our method is en- 
tirely self-contained in the sense that we can compute 
the NMR shielding of both metallic compounds of inter- 
est and the NMR shielding of reference compounds. As 
such, the resulting NMR shifts are directly comparable 
to experimental results. 

The paper is organized as follows. First, we go over 
the physics involved in computing NMR shifts. Secondly, 
we briefly review the so-called "smearing technique^S" 
which allows an accurate and efficient treatment of the 
Fermi surface. We then detail the computation of the 
orbital shift in sec. llVi and of the Knight shift in sec. |Vl 
In sec. I VII we discuss practical issues dealing with the 
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actual implementation of the method. The next section 
is devoted to the study of limit-systems and numerical 
tests. Finally, the last section presents results obtained 
on simple metals. 



II. NMR SHIFTS IN METALS 



field which can be obtained form classical magneto- 
statics 



B.(r') = y J dr <5(r' - r)m(r) 

3m(r) ■ (r' — r) 



dr 



m(r) 



(5) 



A uniform external magnetic field B applied to a 
metallic material generates two different electronic be- 
haviors: (i) a so-called orbital response where electrons 
react to the field as moving charges, (ii) a spin response 
where electrons react as spinning charges. 

In the following and throughout the paper, we use the 
symmetric gauge A(r) = — B A r. with A the gauge, B 
the magnetic field, and r the position in real-space. 

The applied magnetic field induces an orbital current 
Joir'). It can be obtained as the expectation value of the 
current operator J(r'), 



J(r')-J''(r')-^|r')(r'| 
J^(r')--i(p|r')(r'| + |r')(r'|p) 



(1) 
(2) 



c is the speed of light. The first term on the right-hand- 
side of Eq. [l]is the paramagnetic current operator. The 
second term is the diamagnetic current operator as ex- 
pressed within the symmetric gauge. Note that at zero 
field, the expectation value of J(r') is null. 

The orbital current induces in turn an inhomogeneous 
field Bi(r'), which can be obtained from classical mag- 
netostatics, 



B„(r') = i/drjJr)A^ 



r — r 



(3) 



We will describe our approach to the calculation of an 
all-electron induced orbital current using pseudopo- 
tentials in section IIVI The method is an extension to 
metals of the scheme proposed in Ref. |3- 

The spin response results from a spin polarization 
of the electronic cloud by the external magnetic field. 
To compute the resulting net electronic magnetization 
m(r'), we make a co-linearity hypothesis, whereby m(r') 
is supposed parallel to the applied magnetic field. This 
hypothesis is used routinely in hyperfine parameter and 
Knight shift calculations2i2ii^. Hence, m(r') can be ob- 
tained as. 



m(r') = ^h(r')-Pi(0]|, 



(4) 



where p-f(r') and pj,(r') are the up and down spin den- 
sities. The electronic magnetization induces a magnetic 



Bs(r') is composed of two terms: (i) an on-site term, 
called the Fermi contact (first term in Eq. [5]) representing 
the dipole- field at r', (ii) a long-distance dipolar term re- 
sulting from the full magnetic dipole distribution m(r'). 

A method to compute the electronic magnetization 
m(r) to first order in B is given in section Ivl 

For field strengths in the range typical to NMR, the or- 
bital and spin responses can be computed separately and 
to first order in B. The resulting linear relationships be- 
tween the induced first-order fields B^^''(r') and B^^''(r'), 
and the external magnetic field B define the orbital and 
spin shielding tensors, respectively (To(r') and CTs(r'). 

BW(r') = -ao(r')-B ; B« (r') = -a.(r') • B (6) 

The isotropic NMR shielding (t(R) of the nucleus at posi- 
tion R is given by the trace cr(R) = Tr[(To(R) + CTs(R)]/3. 
The isotropic NMR shift S, i. e. the experimental ob- 
servable, is obtained with respect to the isotropic 
shielding dre/ of a so-called zero-shift compound, with 
J(R) = ~(a(R)-a,e/). 



A. Pseudo-potential System 

Within a pseudo-potential system, one must define the 
Hamiltonian and operators with care. Following the pro- 
jector augmented wave methodic (PAW), and the gauge 
including projector augmented wave method^ (GIPAW), 
the spin-Hamiltonian of a system with a homogeneous 
magnetic field becomes 



Ha — 



B Ar 

2c 



R 



sgn(a)^. (7) 



a indicates the spin-channel, p is the kinetic energy op- 
erator, V^^f the magnetic-field-dependent self-consistent 
potential, and V^' the non-local potential at position 
R. sgn((T) returns ±1 depending on the spin channel. 
g^ = 2.0023193 is the gyromagnetic ratio of the free elec- 
tron. The bar above quantities such as Hg- indicates 
pseudo-potential reconstructed operators. The above can 
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be expanded to first order in B as, 



,(0) 



2c 



R 

R / 



HW=sgn(a)B(|+viV(r) 



(8) 
(9) 

(10) 
(11) 



Note that we are interested in systems which are spin- 
degenerate at zero field, hence Ti^'^-' is defined indepen- 
dent of spin and docs not carry a spin index, v^' is a 
reconstruction termi defined as 



(12) 



Square brackets indicate a commutator. L = r A p is 
the orbital momentum operator. In the expansion above 
first order terms are separated into a spin dependent per- 
tubation Tii^ and an orbital dependent term T-L^o^ . The 
former is given within the colinear hypothesis discussed 
in the previous section. 7^1^ is the only spin-dependent 
term in the expansion to first order in B of 7^^- vf;!! 
is the linear part of the self-consistent potential with re- 
spect to B. It is obtained from the functional derivative 
of Vg^f with respect to the first-order electronic magne- 
tization to'^' at zero field. 



2 am(r) 



Hr), (13) 



In order to obtain all-electron NMR shifts, we should 
also reconstruct the current operator and the electronic 
magnetization. The former can be expressed to first order 
as m Ref . i. 



J(r') = J^°^(r')+J'''(r') + 0(S2 



(1), 



with 



and 



J^'^(r') 



j("'(r') = JP(r') + 5^AJ^(r'), 

R 



(14) 



(15) 



B Ar' 



r')(r'| 



E 

R 



2c 

AJi(r') + ^[RAR.r,AJ^(r')] 



(16) 



The paramagnetic reconstruction operator AJ^(r') and 
diamagnetic reconstruction operator AJ^(r') are defined 
as follows. 



AJ^(r') = E IPi^-)[(*R,"|J"(r')l*R. 



(4>R,n|J^(r')l$R,,n)](pR,. 



(17) 



and 

AJ^(r') = _ BA(r'-R) ^ [($R,„|r')(r'|$K,„> 
Zc ^ — ^ 

n.m 

- (^R,n|r')(r'|'i>R,™)](pR,„|. (18) 

The projector functions |pr,,i) are defined in Ref. and 
satisfy (pR,„|0R',m) = (SR,R'(^„,m, where ||0r,„)| is a 
set of pseudo partial-wavefunctions corresponding to the 
all-electron partial wavefunctions {|0R,n)}- 

The electronic magnetization operator ]VI(r') is recon- 
structed using PAW—, 



m{v')^Y.^gn{a)^l^v'){v'\ 

+ Y.\P^M'^^,n\v'){v'\'^TL,. 



I • (19) 

To linear order in _B, the spin and orbital response are not 
coupled. Hence M(r') can be reconstructed using PAW 
only, rather than the gauge including method GIPAW. 

III. METALLIC SYSTEM 

In order to treat the Fermi surface accurately and effi- 
ciently, we follow Ref. [lO and introduce a fictitious tem- 
perature T into the electronic system. 

Let l^j-""*) and ef^ be the eigenvectors and eigenvalues 
of the Hamiltonian 7^*^°^ defined in Eq. [9l Let f{x) be a 
smooth step- function. The occupation fp^i of energy level 

i is defined as fF,i ~ /( — " t ^ )' "^^ere tp is the Fermi 
energy. The latter is recovered from the conservation of 
the number 27V of electrons in the system, ^ ■ fp.i = 2N, 
with i running over all eigenstates. 

It was shown by de Gironcoli^ that the first order 
expectation value o'-'^' of an operator O = C''^-' + 0'--^\ 
with ((^)) indicating the zero (first) order pertubation 
expansion, can be recovered as 

-f2E/F,.(^fV'l^r^) 



2E^^( 



ji) ,(0) m 



T 



)(§(0)|O(0)|^(0)). (20) 



5R is the real value. The sum over i runs over all states. 
1-6'^^ is some pertubation (it will be either the orbital or 
spin pertubation of Eqs. [TO] and [TT1) . The last term ac- 
counts for variations of the Fermi energy to first order 
The linear variation of the Fermi energy can be 
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recovered from the conservation of the number of elec- 
trons, J2iT~\e^F - £f^)(5(- '''^°'~''°' ) = 0. In this work 
we will always have e^' — 0. Function S{x) is defined 
as the derivative of f{x), S{x) = — df(x)/dx. The Green 
fimctions Q{e) is defined as 

^(^o=E%f-%i^r)(^ri (21) 

The sum over j runs over all states. For i = j, the 

limit ifp^, - - ef^) '3 (5(-^2^) is taken. 

Expression 1201 contains a factor two for spin. 

IV. NMR ORBITAL SHIFTS 

The method presented in this section is an extension 
to metals of the scheme proposed in Ref. to compute 
NMR shifts in insulators. The Fermi surface is modeled 
using the smearing scheme of Ref. H^. For the sake of 
simplicity, the proof is given for an all-electron system 
(i. c. with ^ 0). 

We first compute the induced-current to first order for 
a finite system. The result is rc-cxprcssed in a form suit- 
able for extended systems using the sum-rule of appendix 
lAl This expression is then specialized to the case of pe- 
riodic systems. Finally, we give the expression of the 
orbital current for a pseudo-system. 

A. Finite systems 

By setting = 0, the Hamiltonian of an all-electron 
system is recovered from Eq. [8l 

H^n^°^ +H'^^^ +0{B^), (22) 
HW = ip2 + v.,/(r) (23) 

7^(1) = ■ B (24) 

We note |^,) the all-electron wavefunctions. The cur- 
rent operator for an all-electron system is given in Eq. [T] 
Using the linear response Eq. [20l the expectation value 
j'"'^-'(r') can be recovered as, 

j(i)(r') = 2^5R{(*(°)|J^(r') ^(.f)) H^^l^r^)} 

i 

_^,(0)(r'). (25) 
c 

In the above equation, we have used the assumption that 
there is no linear order variation of the Fermi energy, 
e'p' ~ 0. Indeed, in a non-degenerate system, the linear 
order variation of the eigenvalues are t^P = {^f^\j^'L A 
B|5'^°'') for a given field B. Since the zero order system 



is invariant upon time reversal, the wave-functions |^,- } 
can be chosen real. Hence, we have el^"* = 0. It follows 
from the condition on e^' given in section Hill that = 
0. 

Eq.[2S]is valid for a finite system only, indeed, for r' i— > 
oo, BAr'p(r') diverges in an extended system. There is a 
similar divergence in the other term of Eq. 1251 such that 
the orbital current itself is finite. Yet, from a numerical 
point of view. Eg . [25] cannot be used to compute j^^-* (r'). 

B. Extended System 

Following Ref. 0, Eq. [25] can be reexpressed using a 
generalized /-sum rule (given in appendix [X|) into a more 
practical expression for an extended system. We have. 



^/^,,(vl/f |-[BAr'.r,J^(r)]|vI/(°)), (26) 

i 

where is an odd operator and BAr'-r an even operator. 
Using the sum rule, Eq. [25]can be rewritten as, 

j(i)(r') = iE^{(*i''^|J^(r') 

c;(ef)) (r-r')Ap|*f')| (27) 

Since position quantities now enter as differences, it fol- 
lows that the above expressions is invariant upon trans- 
lation of the system. Furthermore, the Green function at 
finite temperature is short-ranged. It follows that contri- 
butions to the orbital current vanish for large values of 
(r - r') in Eq.[22] 



C. Periodic System 

At this point, we have a formalism adequate for obtain- 
ing the current response in extended metallic systems. Of 
those, only translationally-invariant periodic systems are 
computationally feasible. Hence, we now introduce these 
translational symmetries explicitly into the equations for 
the current response. We write 1^''^'') = ^''^ '^W^k) 
electronic Bloch states of crystal momentum k. is 

the corresponding eigenvalue. (r|M^|^^) is a normalized 
cell-periodic function. In the spirit of Ref. we trans- 
form the real-space dependence (r — r') into a reciprocal 
space dependence by introducing the limit, 

(r - r') lim — V L^''^--^"-'"') - e-""""-^"-"")] , 

OL—x,y,z 

(28) 
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where Ua=x,y,z is real-space basis. This transformation is 
subject to the condition |r — r'| < C (C a vector) which 
is verified since contributions to the orbital current in 
Eq. [57] vanish for large values of (r — r'). The orbital 
current is then recovered as a numerical derivative, 



lim-[S(r',g)-S(r',-q)] 
q^o 2q 



(29) 



where, 



a=x,y,z i,k k 



"ik |Jk,k+gu„('' ) 



ak+9u„ (e,k) B A u„ • (p + k)\u\l^) j . (30) 

Nk is the number of k-points in the discrete integration of 
the Brillouin zone. Wc have introduced the k-dcpendcnt 
Green function Gk{i), 



5.ie) = E (V " (31) 

and the k-dependent paramagnetic current operator 
IP 

"^k,k" 

Jk.k'(r') = -I (P + k) |r')(r'| - i|r')(r'| (p + k') (32) 

Eq. [29] allows us to compute the orbital current of an 
all-electron system. In practice, it is more efficient to 
use pseudo-potentials when expanding the density on a 
plane-wave basis set. We now give a general expression 
for the orbital current in periodic pseudo-potential sys- 
tems using the GIPAW reconstruction scheme of Ref. |j. 



D. Periodic pseudo-potential system 

The orbital current can be obtained from apseudo- 
system using Eq. [5] and Eq. [T3] Following Ref. |j as well 
as the steps given above, one can find an expression for 
the orbital current suited to a periodic pseudo-system. 

We find that the current is composed of three compo- 
nents: (i) the bare current j[^'^g(r'), (ii) the paramagnetic 
augmentation current j^p(r'), (iii) the diamagnctic aug- 
mentation current j Ad 



j^'^(r')=£L(r')+jk'^(r')+j«(r') 



(33) 



The diamagnetic augmentation current is simply the 
expectation value of the operator given in Eq. I18i 

ji'kr') = 2E(*!k^|AJi(r')|*l^^) (34) 



Note that the projectors |pn,R) rnake AJ^(r') short- 
ranged. Furthermore, since positions quantities enter as 

differences, jAd(''') i^ translationally invariant. 

The paramagnetic augmentation and bare currents are 
obtained as numerical differences, 

jiale(r') = lim ^ [Sba.e(r', q) ~ Skareir', ~q)] , (35) 



ji^J(r') = lim ^ [SAp(r', q) - SAp(r', -q)] . (36) 
The two newly introduced functions are defined as, 

S,_(r',g) = ^ J2 J:^\'-{u^'\Km+,^A^') 

0!. — x,y,z i,k: L 



(5k+9UQ(eik) B A Uq • Vk+qu„,k|Mik ) f"' (37) 



and, 



a — x,y,z i.k V 

5k+quQ(ejk) B A Uc ■ Vk+gu„,k|Mik )|- (38) 

\u[^) is the cell-periodic function such that l^'^^ ) = 
gjkAr|y(o)^^ The Green function ^k(e) is redefined using 
the pseudo-eigenstatcs. 



fpjk - / - 



(0) 



(39) 



A k-dependent non-local pseudopotential V^'^., is also 

defined, which acts on k Bloch states on the left and k' 
states on the right, 



(40) 



T n, .771 



The periodic projectors „) are obtained from the real- 
space projectors |pL+T-,n) as 



where the sum runs over the lattice vectors L. Cell- 
internal atomic-coordinates are noted with r. The ve- 
locity operator is also redefined as, 

Vk,k'=P + k' + i[r,V^,'k']- (42) 

Finally, a k-dcpendcnt paramagnetic current operator 
k' affiliate pseudo-operator AJ^ ^ k' ™" 

troduced. 

Jk,k'(r') = (P + k) |r')(r'| - lW){r'\ (p + k) (43) 
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n,m 

-(^L+..„|J^(r')|0L+.,™)](p'?,„| (44) 

The orbital shielding is then obtained from Eq. [3] and 
from its definition Bq — — CTo • B. 



shifts in metals. It was implemented in a parallel plane- 
wave pseudopotential electronic structure code. We now 
outline the features specific to the NMR method. We 
shall first discuss the application of the Green function, 
common to both orbital and Knight shift computations, 
and then turn to the specifics of each type of response. 



KNIGHT SHIFT 



A. Linear response 



We now turn to the Knight shift, which results 
from the electrons interacting with the field as spinning 
charges. More specifically, the magnetic field induces a 
net electronic-spin which then interacts with the mag- 
netic nuclear dipole through the Hyperfine interaction 
(Eq. [5]). The Knight shift measures this interaction. 

The Hamiltonian to first order is given up to first order 
be Eq. Eq. H and Eq. \n\ 



1 



V, 



(0) 
scf 



H«=sgn(.)B(|+viV 



(45) 
(46) 



The linear order wavefunctions l^",-^^) and eigenvalues e 
are anti-symmetric with respect to field direction, i. e. 



when B is mapped onto B 



-B, we expect e. 



(1) 



.(1) 



and \^ 



^ts)^ where cr is the spin opposite 



to cr. It follows then that \^)^') = -l^^^) and e,^| 
^^ti ■ From this last condition, it follows that there is 

no variation of the Fermi energy to first order, e^P = 0. 

For simplicity, the following is obtained directly for the 
pseudo-system. Indeed, the reconstruction of the con- 
stant part of Ti!'a^ is zero. Furthermore, we neglect the 
polarization of the core electrons by the valence spin- 
density. In practice, this is equivalent to neglecting the 
PAW reconstruction of the self-consistent pertubation. 

Exploiting the spin anti-symmetry described above, 
the electronic magnetization to fist order in B can be 
obtained as, 



.(1) 



m(i)(r') = 2^3?i (^-f ^ |M(r') 



(47) 

with the quantities defined previously. 

Once the electronic magnetization is obtained, the 
Knight shift can be computed from Eq. [5] and Eq. [6] 



/(I) 



(1) 



|(0)n 



VI. PRACTICAL IMPLEMENTATION 

The goal of the method presented above is to provide a 
practical and quantitative approach to computing NMR 



we are interested in computing first-order quantities 
(see Eqs.[35l[36l and[47| such as. 



.(1) 



(48) 



where O is an operator and Ti*-^-* some pertubation. the 
green function is expressed as in Eq. [2TJ Both the sum 
over i above, and that over j in Q{e) range over all states, 
such an expression cannot be calculated directly, it was 
shown by de Gironcoli in ref. [13 that o'^' can be com- 
puted via an alternate first-order wavefunction 



(49) 



such that the sum over i runs only over partially occupied 
states. I^^f^) can be also computed without reference 
to empty states. 



7^(0) + Q 



JO) 



(0): 
3 I 



nT 



Q £ , r I fp.i If,] 

Pi,3 = jF,t9i,] + jF.j9j.i + Oij—) i0)9],: 



(50) 



g{x) is a symmetric function such that g{x) -\- g{—x) = 1. 

We define gij = g( ' ^ ' ). Partially occupied wave- 
functions are defined such that ef^ < — nT < 
{ai 7^ 0), where n is a suitably large number. We find 
that orbital and spin shieldings are converged for n = 7. 



B. Orbital shifts 

The method presented above differs only slightly from 
the prior method for insulators. We will address only 
these differences and defer the interested reader to Ref. [j- 

The macroscopic induced field B(,^'(G = 0), where G 
is a vector of reciprocal space, is not a bulk property. 
Indeed it results from the surface current in the sample, 
and hence depends on the shape of the sample. Fol- 
lowing Ref. 0, we compute it through the so-called bare 
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macroscopic susceptibility Xbare, consistent with the on- 
site approximations for the reconstruction current, 

BW(G = 0) = ^4^Xfcare-B, (51) 

Xbare is the Contribution to the macroscopic susceptibil- 
ity from the bare current We adapt the ansatz of 
Ref. to the case of metallic compounds, 



Xbare = Hm ^ 



(52) 



F{q)+F{-q)-2F{Q) 
where Fij = {2 — 6ij)Qij{q). i and j are Cartesian indices 

a — x.y,z i.k 



{u\^\Ua A (p + k) C/k+qu„(e-k ) Ua A Vk+qu„,k|u. 



{0)\ 



(53) 



When interested specifically in the susceptibility Xo, 
we use another ansatz from Ref. 0, with 



Xo = lim ^ F*°*(g) -t- F'°'{^q) ~ 2F*°*(0) 
<?i-»o q^ L 



(54) 



a— a:.y,z i.k 



/-(0)l 

{ul^^ \Ua A Vk,k+gu<:, 

ek+,ju„(e*k ) u,, A Vk+,u„,k|u'k > [ (55) 



and = (2 - J,,)Q*f (g) 

At zero temperature, Xo and Xbare above and the 
corresponding quantities of Ref. are equivalent. 



C. Knight shift 



The variation of the self-consistent potential vf (r) is 
evaluated using a simple self-consistent loop over the cal- 
culation of the first-order wave-functions. In other words, 
the spin density is recomputed at each step and V^^|(r) 
updated. In the case of local density approximations, 

'^^lc/(i") is simply, 

vLH(r) = ^^P«(r), 



'f'^^ (56) 
(vL(r)-Vi,(r))pW(r), 



dpsir 



where Vj^ and Vj-^ are the exchange-correlation potential 
of the up and down spin channels, respectively, computed 
from the ground state densities. pi^\v) = p'"^^\r)—p^^'' (r) 



is the first order spin density at r. These derivatives are 
evaluated numerically for each point of the real space 
mesh. The self-consistent Hartree potential is not spin 
dependent, and hence it is not modified by variations of 
the spin density. 

A pertubation using generalized gradient approxima- 
tions can be implemented in much the same way. 

We find that convergence with respect to the number of 
iterations over V^^j-(r) can be achieved efficiently without 
mixing. 



VII. NUMERICAL TESTS 

A. Interacting Homogeneous gas 

NMR shifts require the computation of the macro- 
scopic susceptibility in order to account for the diamag- 
netic shielding resulting from surface currents. We will 
now test these calculations against available analytical 
results for the homogeneous electron gas. The orbital 
(Xo) E^nd spin (xs) susceptibilities per unit volume of this 
model system are given by the formulae^: 



1 



Xo 
Xs 



12c2 

9e 



1 



(57) 



.9(eF) 



3 N 



1/3 



N is the number of electrons in the system, exc is the 
exchange-correlation energy per unit volume as given by 
PBE^i, and giep) is the density of states at the Fermi 
energy. The derivative of e^c is evaluated numerically. 

The fractional factor in Xs results from the exchange- 
correlation. More specifically, the magnetic field in- 
duces a polarization of the electrons at the Fermi en- 
ergy, which then propagates to lower lying levels through 
exchange-correlation interactions. Indeed, for a non- 
interacting homogeneous gas, the spin susceptibility re- 
duces to Xs = 3e/ (4c2)5(ei?), i. e. it is simply proportional 
to the available degrees of freedom at the Fermi surface. 
This propagation effect is rendered computationally by 
the self-consistency of Eq. HTl 

To simulate a homogeneous gas within a pseudo- 
potential code, we construct a pseudopotential with zero 
potential and zero atomic charge. A temperature of 
0.4 eV is introduced into the system. We use an fee unit 
cell with a cell-parameter of 3.6lA. The Brillouin zone is 
sampled with a 60x60x60 Monkhorst-Pack grid. Different 
electronic densities are obtained by varying the number 
of electrons in the cell. 

Results are given for a range of densities (parameter- 
ized by rs/uo = (3/477/5)"^, where ag is Bohr constant) 
in Fig. [TJ X dots represent the response computed with 
our approach, and solid lines are analytical results. The 
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PBB^ exchange-correlation functional is used. Results 
agree to within numerical noise. 
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we use a spherical sample when accounting for surface 
currents. We use the PBE^i exchange-correlation func- 
tional. 

Aluminum and copper are cubic face centered metals 
with a ~ 4.05 A and a ~ 3.61 A, respectively. Lithium is 
body centered with a = 3.49 A. We use experimental cell 
parameters as given by Ref. [H. 



Metal 


Smearing 


Cutoff 
Knight Orbital 


Knight 


Afc 

Orbital 


Al 


0.15 


15 15 


29820 


62790 


Li 


0.2 


15 15 


8094 


11900 


Cu 


0.2 


75 90 


1300 


5740 



TABLE I: Computational Details. Smearings are given in 
eV. Plane wave cutoffs for both Knight and orbital shift cal- 
culations are given in Rydberg. Nk stands for the number 
of independent k-points in the irreducible wedge of the Bril- 
louin zone. The Brillouin zone is represented with a discrete 
Monkhorst-Pack grid^*^. The Marzari-Vanderbilt smearing 
functions^ is used. 



B. Zero shift compounds 



FIG. 1: Spin and Orbital susceptibility per unit volume of 
the interacting-electron homogeneous gas respect to Vs/ao = 
(S/iivp)^. Solid lines represent analytical results while dots 
represent results computed with this code at 0.4 eV smear- 
ing. The susceptibilities are dimensionless. The exchange- 
correlation is modeled with the PBE functional. 



VIII. SIMPLE METALS 

The object of the present work is to build a quan- 
titative method for computing NMR shifts in metallic 
compounds. As such, we now study three simple metals: 
bulk aluminum, bulk lithium, and bulk copper. 

Experimentally, NMR shifts are obtained with respect 
to the response of so-called zero-shift compounds. We 
will first study this aspect of the problem, and compute 
the shielding of these compounds. We will then give the 
computational details for each metal, and finally examine 
the NMR shifts and macroscopic magnetic susceptibili- 
ties. 



A. Computational Details 

Computational details are reported in tables [J For 
all calculations, we use the Marzari-Vanderbilt smear- 
ing function^^ and TrouUicr-Martins^ norm-conserving 
pscudopotentials. Following experimental conventions. 



Experimental NMR shifts arc obtained as 

U • (5 • U = — U • (<T — f?re/) • U, 



(58) 



where u is the direction in which the external magnetic 
field is applied, (t is the the shielding of the compound, 
and (T, e/ is the shielding of the zero-shift compound. 

Rather than evaluating (T^e/ directly, we will compute 
the shielding of some compound for which the NMR shift 
is well known experimentally, and then deduce CT^e/ from 
Eq.EHl 



Atom 


"zero-shift compound" 


compound 










type 




^exp 




Al 


AICI3 in heavy water 


AIPO4 


519 


45 


564 


Li 


aqueous LiCl 


LiaO 


86 


10 


96 


Cu 


CuBr powder 


CuBr 


424 


0.0 


424 



TABLE II: Reference shifts a^'^-^ . The "compound" columns 
gives the solid which is used to obtain the reference shift, 
its calculated isotropic shielding cr"*, and its experimental 
isotropic shift S^xp- The reference shift is obtained using the 
relationship Sexp = — (cr"' — a^^-^). Shieldings are converged 
to better than a ppm. Shieldings and shifts are given in ppm. 



The reference shifts for each element Al, Li, and Cu are 
given in Tab. [Ill Note that all references are computed 
on insulators, and hence that the shieldings result only 
from the orbital response. The latter are computed using 
the method for insulators described in Ref. \ 
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C. Behavior with respect to smearing 

The computation of NMR shifts requires a very fine 
description of the Fermi surface. Hence, one must take 
care that the computed shifts are indeed converged with 
respect to smearing. Figures [5] and [3] report the con- 
vergence behavior with respect to smearing of, respec- 
tively, the spin macroscopic spin-susccptibihty, and of 
the macroscopic orbital-susccptibihtics for Aluminum, 
Lithium, and Copper. Figures [3] and O report the behav- 
ior of the Knight shift and of the orbital shift, excluding 
the contribution of the macroscopic susceptibility. We 
find that the orbital susceptibility is the hardest to con- 
verge. This is coherent with the fact that as a second 
order derivative of the total energy, it depends on very 
fine details of the Fermi surface. On the other hand, 
the spin susceptibility is obtained as the average over the 
unit cell of the spin density. As such, it is comparatively 
insensitive to details of the Fermi surface, and converges 
much faster with respect to the smearing parameter. A 
similar hierarchy is obtained for the convergence behav- 
ior of the Knight and orbital shifts (not including their 
respective susceptibility). It should be noted that in the 
examples provided here, the Knight shift is by far the 
largest component of the total NMR shifts. Overall, we 
expect the total NMR shielding to be converged to better 
than 4% with respect to smearing and k-point density. 

On the other hand, convergence of the magnetic sus- 
ceptibility can prove quite arduous. For instance, the 
orbital susceptibility of Aluminum varies from —0.3 
to -f 5.6 10~^ cm'^ mol^^ within the temperature range 
0.3 eV to 0.1 eV. Aluminum presents the slowest conver- 
gence of the three metals studied in this work. 

D. Results and Discussion 

1. Macroscopic Magnetic Susceptibility 

The computed magnetic susceptibility are referenced 
in Tab. IIIII Overall, agreement is very good. It contains 
a diamagnetic contribution from the core electrons. This 
contribution is constant within the frozen core approxi- 
mation and is computed once and for all from an atomic 
code for each pseudo-potential. Tab. IIVI compares the 
spin and orbital susceptibilities of each metal to an elec- 
tronic gas of corresponding mean density. 

When examining the band structure of Aluminum, one 
finds that it is quite similar to that of an homogeneous 
gas of equivalent density. As a result, the non- interacting 
spin susceptibility and the Stoncr factor of these two sys- 
tems are remarkably close. This indicates that not only 
are their density of states at the Fermi level similar, but 
also the Pauli-mediated behavior of the electrons with 
respect to a pertubation of the spin population. On the 
other hand, the orbital susceptibility of these two systems 
are quite different (note however that for Aluminum, we 
did not achieve good convergence of this quantity with 




Smearing (eV) 



FIG. 2: Convergence with respect to smearing (a) of the spin 
susceptibility of Aluminum (green triangles). Lithium (red 
squares), and Copper (blue +). For comparison purposes, the 
spin susceptibility of each metal is normalized to its value at 
the lowest achieved smearing. The a;-axis represents smearing 
in eV. 



respect to smearing). Indeed, in an ideal gas, the contri- 
bution of lower lying electrons cancels-out exactly. Thus, 
only electrons at the Fermi level contribute to the orbital 
susceptibility. This is usually not true in more complex 
systems. Even small differences between the band struc- 
tures of Aluminum and the homogeneous gas will result 
in appreciably different orbital susceptibilities. 

Lithium presents a case very different from the one 
above. Its non-interacting spin susceptibility is much 
larger than that of the homogeneous gas. As a result, 
the large polarization at the Fermi level yields a large 
polarization of the lower-lying electronic wavefunctions. 
The Stoner factor of Lithium is much larger than that 
of the homogeneous electron gas. Interestingly, Lithium 
presents very little orbital susceptibility. 

Copper presents a different picture still. Indeed, it 
has a rather low density of states at the Fermi level 
compared to the homogeneous gas. As a result, both 
non-interacting and interacting spin-susccptibilitics are 
small. On the other hand, the large number of lower ly- 
ing electrons, including d electrons, yield an appreciable 
diamagnetic orbital susceptibility. As such, of the three 
metals studied here, it is the only one with a diamagnetic 
susceptibility. It is worthwhile to note that only the or- 
bital susceptibility can explain such a behavior, and that 
hence a complete understanding of the susceptibility of 
Copper requires the computation of both spin and orbital 
contributions. 
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0.2 0.4 0.6 0.8 1.0 
Smearing (eV) 

FIG. 3: Convergence with respect to smearing (a) of the or- 
bital susceptibility of Aluminum (green triangles), Lithium 
(red squares), and Copper (blue +). For comparison pur- 
poses, the orbital susceptibility of each metal is normalized 
to its value at the lowest achieved smearing. The a;-axis rep- 
resents smearing in eV. 



Metal 


Xs 




Xo 




\core 


x'' 


Exp. 


Al 


17.7 ± 


0.2 


1.9 ± 


5 


-3.0 


16.6 


16.5 [28] 


Li 


28.4 ± 


0.5 


0.7 ± 


1 


-0.7 


28.4 


24.5±0.3 [29] 


Cu 


10.8 ± 


0.2 


-13.1 ± 


1 


-4.5 


-6.8 


-5.3 [30] 



TABLE III: Isotropic magnetic macroscopic susceptibility (in 
10^" cm'^mol"^, moles of atoms). The susceptibility con- 
tains three components: (i) the orbital susceptibility (xo), 
(ii) the spin susceptibility (xa), and (iii) the diamagnetic sus- 
ceptibility of the core electrons (Xcore). As shown in Fig. [S] 
we were unable to converge the orbital susceptibility of Alu- 
minum. 



2. NMR shifts 

The computed isotropic NMR shifts arc reported in 
Tab. |Vl Tab. I VII also report as/cr^g, a quantity akin to 
the Stoner factor of the susceptibility, where is the 
Knight shift computed including self-consistency, and cr^ 
the Knight shift computed without self-consistency. 

The NMR shift of Aluminum results predominantly 
from the Knight shift. It is worthwhile to note that the 
orbital and Knight shielding tensors are of similar magni- 
tude, — 548ppm and 1862 ppm respectively. Yet, whereas 
the Knight contribution enters into the NMR shift as a 
whole, the orbital part enters as a variation of the abso- 
lute orbital shielding tensor between pure Al and ionic Al 
(which presents no Knight shift), yielding a much smaller 
contribution. Previous theoretical calculations^ predict 




Smearing (eV) 



FIG. 4: Convergence with respect to smearing (a) of the 
Knight shift (not including the spin susceptibility) of Alu- 
minum (green triangles). Lithium (red squares), and Copper 
(blue -I-). For comparison purposes, the Knight shift of each 
metal is normalized to its value at the lowest achieved smear- 
ing. The X-axis represents smearing in eV. 



System 




Stoner 


Xs 


Xo 


Al 


13.2 


1.34 


17.7 


1.9 


gas 


12.5 


1.31 


16.4 


-4.2 


Li 


15.5 


1.83 


28.4 


0.7 


gas 


10.2 


1.48 


15.1 


-3.4 


Cu 


9.5 


1.14 


10.8 


-13.1 


gas 


15.3 


1.18 


18.1 


-5.1 



TABLE IV: Isotropic magnetic macroscopic susceptibilities 
(in 10~^ cm^mol"^, moles of atoms). Xs is the non- 
interacting spin susceptibility, Xs the interacting spin suscep- 
tibility, and Xo the orbital susceptibility. For comparison, the 
susceptibilities of a homogeneous gas of the same mean den- 
sity as the system is given. 



a Knight shielding of CTs = 1707 ppm. Although, the au- 
thors of Ref. [m do not compute NMR shifts comparable 
to experiment, in the sense that they do not reference 
their results to a computed zero-shift compound, their 
result is close to experimental value because of the pre- 
dominance of the Knight shift. As will be the case for the 
other metals studied here, the ratio (Js/cr^ and the Stoner 
factor are quite close in value. Indeed, both quantities 
represent the same physical phenomena, namely the in- 
terplay between the Kohn-Sham potential of the valence 
electrons and the spin polarization at the Fermi level. 

Again, the orbital shift of Lithium is by far smaller 
than its Knight shift. As mentioned previously, the lower 
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0.2 0.4 0.6 0.8 1.0 
Smearing (eV) 

FIG. 5: Convergence with respect to smearing (a) of the or- 
bital shift (not including the orbital susceptibility) of Alu- 
minum (green triangles), Lithium (red squares), and Copper 
(blue -I-). For comparison purposes, the orbital shift of each 
metal is normalized to its value at the lowest achieved smear- 
ing. The a;-axis represents smearing in eV. 




r H p N r p 



FIG. 6: Band structure of lithium. The width of the line is 
representative of its contribution to the isotropic Knight shift. 
The Fermi energy is set to zero on the energy scale. Compu- 
tations were done with a smearing of 0.2 eV, for which the ^ 
dependence is obvious at the Fermi energy. The "divergence" 
in T disappears with the Brillouin zone integration. 



lying levels are heavily polarized by electrons at the Fermi 
surface. The authors of Ref. [33 estimated the Knight 
shift of Lithium including core-polarization. Even in this 
case, where from Fig. [6] one would expect a rather high 
polarization, the contribution is only of the order of 5% 
of the whole (250 ppm). More recently Mishra et al esti- 
mate a Knight shift of 301.9 ppm^. Overall, our calcula- 
tion agrees very well with experimental values. The ratio 
(Js/<^s is relatively smaller than the Stoner factor. One 
should note that the latter is a ratio of the average spin- 



polarization over the whole unit cell, whereas the former 
is the ratio over the spin-polarization at a single point 
of unit cell, namely the position of the Lithium nucleus. 
The discrepancy between the two quantities implies sim- 
ply that the effect of the spin polarization is smaller at 
the nucleus than on average across the cell. 

Of the three metals studied here. Copper is the only 
one which presents an appreciable orbital contribution 
to the NMR shift. It is probably a result of the filled ri- 
bands. Nonetheless, as large as the orbital contribution 
may be, the Knight shift is larger still. Interestingly, 
the computed absolute orbital shielding tensor (includ- 
ing both valence and core contributions) is rather small 
(26 ppm) . It would seem that a substantial paramag- 
netic contribution from the valence electrons cancels out 
the substantial diamagnetic contribution from the core 
electrons (computed to be 2171 ppm). In other words, 
whereas in Li and Al, the reference compound and the 
metals had similar orbital shielding tensor, the orbital 
behavior of metallic Cu is very different from that of 
Copper-Bromide. As was the case for the magnetic spin 
susceptibility, the spin-polarization at the Fermi level has 
little effect on the lower lying levels, resulting in a rela- 
tively small o's/o's ratio and Stoner factor. 



Metal 


Us 




0"o — CTref 


s 


Exp. 


Al 


-1858 ± 


70 


-16 ± 8 


1874 


1640 [33|, 


Li 


-266 ± 


5 


-15 ± 1 


281 


260 [28| 


Cu 


-2336 ± 


20 


-450 ± 10 


2786 


2380 [28], 



TABLE V: Isotropic NMR shifts of a few simple metals. For 
comparison, the orbital shielding with respect to the reference 
and the Knight shifts are given as well. The isotropic NMR 
shifts are given by the relationship i5 = —{ao + cts — (Jre})- 
Estimates of the convergence with respect to temperature and 
Brillouin zone sampling are given in the first two columns. 



Metal 






<Js 


(To — CTre/ 


5 


Exp. 


Al 


-1330 


1.40 


-1858 


-16 


1874 


1640 [33], 


Li 


-157 


1.69 


-266 


-15 


281 


260 [28] 


Cu 


-2121 


1.10 


-2336 


-604 


2940 


2380 [28], 



TABLE VI: Isotropic NMR shifts of a few simple metals, a], is 
the non-interactmg Knight shift computed without the self- 
consistent part of the pertubation. The ratio (Js/cr1 is the 
Knight shift equivalent of Stoner factor of the spin suscepti- 
bility. Unsurprisingly, this ratio is quite close to the Stoner 
factor. Indeed both are a measure of the interplay between 
the spin-polarization at the Fermi level and lower-lying va- 
lence electrons. 
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IX. CONCLUSIONS 

We have presented a unified method for computing 
NMR shifts in metals. Our approach yields shifts which 
are directly comparable to experimental data, in the 
sense that both orbital and Knight shifts are computed. 
It was implemented within a pseudo-potential, plane- 
wave density functional theory code. All-electron quan- 
tities were recovered using the PAW approach. Gauge 
invariance was enforced with GIPAW. We compared re- 
sults given by our approach to known analytical solu- 
tions for the homogeneous gas. Finally we successfully 
computed the NMR shifts of simple metals, with good 
comparison to experimental results. In conclusion, we 
have described a method which can accurately recover 
the NMR shifts of real metallic systems, thus allowing a 
better interpretation of NMR data. Next, we expect to 
study semi-metallic systems, such as graphite and nan- 
otubes, for which an accurate description of both orbital 
and Knight shift is of paramount importance. 
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APPENDIX A: THE GENERALIZED /-SUM 
RULE 

Let O and £ be odd and even operators respectively 
on time reversal, i.e. for any real wave- functions j^*) and 
I*'): 

(*|0|*') = -(*'|C|*)> I*') = (Al) 

Let l^'j) be the eigen- wave- functions of the haniilto- 
nian Ti, with eigenvalues e^. Let f{x) be a smearing 
function and cr the smearing. Then the occupation fac- 
tors are defined as fj.i = f{-^-^) (where i,j = F stands 
for the Fermi energy ep, and finally, let 

s^Y. 'fi{{^^\Og{e.)][£,n]\^A (A2) 

c;(e,) = ^/^^l-i^|vI/^.)(vI/^-| (A3) 



where 5R is the real part. Then, using the fact that 
TLY^i) = we arrive at the expression: 



1 



which can be separated into two sums: 



{^,\0\^,){^,\-£\^.,) 



1 



{^.\0\^,){^,\-£\^,) 



Swapping dummy indexes in the second term: 



{^.,\0\^,){^,\-£\^,) 



{^,\0\^,){^.\-£\^,) 



Then, using the parity of O and £: 

1 



{^.,\0\^,){^,\-£\^.) 



[■^,\0\^,){^j\-£\^.,) 



After remarking that X), l = 1: 



{^,\-0£\^,) 



(A4) 



(A5) 



(A6) 



(A7) 



(A8) 



s = 2E./F,«3ft 

i 

Expanding the real value, we arrive at the result: 
-E '^[{^^\Og{e,)\[£,H\\'^^)\ = 

Y,fp,{'^^\][^M\^^) (A9) 



Expression IA9I and equation (A7) in the appendix of 
Ref. 3 differ by the definition of the Green function and 
the range of the sum over states. At zero temperature 
and in insulators, the results arc equivalent. 
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